home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_5 / a5_5.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.3 KB  |  141 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 5.5 (Trigonometric Polynomials).
  9. % Section    5.4,    Fourier Series and Trigonometric Polynomials, Page 311
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % This program finds the trigonometric polynomial
  13.  
  14. % for periodic data in the interval [-pi,pi].
  15.  
  16. % Remark. tpcoeff.m and tp.m are used for Algorithm 5.5
  17.  
  18. pause % Press any key to continue.
  19.  
  20. clc;
  21. %    The computer will place the abscissas in  X
  22. %    The computer will place the ordinates in  Y
  23. % Place the f(x) in the 'string function'  fun
  24. % Place the number of points in  n.
  25. % Place the order of the trig. poly. in  m
  26.  
  27. fun = 'x/2';
  28. n = 12;
  29. m = 5;
  30.  
  31. h = 2*pi/n;
  32. X = -pi:h:pi;
  33. x = X;
  34. Y = eval(fun);
  35.  
  36. [A,B] = tpcoeff(X,Y,m);
  37.  
  38. pause    % Press any key to find the trigonometric polynomial.
  39.  
  40. clc;
  41. hs = 2*pi/150;
  42. Xs = -pi:hs:pi;
  43. Ys = tp(A,B,Xs,m);
  44. a = -3.5;
  45. b = 3.5;
  46. c = -2;
  47. d = 2;
  48. axis([a b c d]);...
  49. plot(X,Y,'or',Xs,Ys,'-g');...
  50. hold on;...
  51. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  52. xlabel('x');...
  53. ylabel('y');...
  54. Mx1 = 'Trigonometric polynomial: y = P';...
  55. Mx2 = [Mx1,num2str(m),'(x)'];...
  56. title(Mx2);...
  57. grid;...
  58. axis;...
  59. hold off;...
  60. shg; pause    % Press any key to continue.
  61.  
  62. clc;
  63. format short e;
  64. Mx1 = 'The trigonometric polynomial has been determined.';
  65. Mx2 = 'The degree is  m = ';
  66. Mx3 = [Mx2,num2str(m)];
  67. Mx4 = 'The coefficients a(n) for cos(nx) are:';
  68. Mx5 = 'The coefficients b(n) for sin(nx) are:';
  69. clc,echo off,diary output,...
  70. disp(Mx1),disp(''),disp(Mx3),disp(''),...
  71. disp(Mx4),disp(''),for i=1:5:m+1,disp(A([i:min(i+4,m+1)]));end,...
  72. disp(Mx5),disp(''),for i=1:5:m+1,disp(B([i:min(i+4,m+1)]));end,...
  73. diary off,echo on
  74.  
  75. pause % Press any key to continue.
  76.  
  77. clc;
  78. format short;
  79. points = [X;Y;tp(A,B,X,m);Y-tp(A,B,X,m)]';
  80. Mx6 = '    x(k)      y(k)      P(x(k))   error';
  81. clc,echo off,diary output,...
  82. disp(''),disp(Mx6),disp(points),diary off,echo on
  83.  
  84.  
  85.  
  86. pause    % Press any key to find another trigonometric polynomial.
  87.  
  88. clc;
  89. % Place the order of the trig. poly. in  m
  90.  
  91. m = 2;
  92.  
  93. [A,B] = tpcoeff(X,Y,m);
  94.  
  95. % Remark. Now the degree is only m = 2.
  96.  
  97. pause % Press any key to continue.
  98.  
  99. hs = 2*pi/150;
  100. Xs = -pi:hs:pi;
  101. Ys = tp(A,B,Xs,m);
  102. a = -3.5;
  103. b = 3.5;
  104. c = -2;
  105. d = 2;
  106. axis([a b c d]);...
  107. plot(X,Y,'or',Xs,Ys,'-g');...
  108. hold on;...
  109. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  110. xlabel('x');...
  111. ylabel('y');...
  112. Mx1 = 'Trigonometric polynomial: y = P';...
  113. Mx2 = [Mx1,num2str(m),'(x)'];...
  114. title(Mx2);...
  115. grid;...
  116. axis;...
  117. hold off;...
  118. shg; pause    % Press any key to continue.
  119.  
  120. clc;
  121. format short e;
  122. Mx1 = 'The trigonometric polynomial has been determined.';
  123. Mx2 = 'The degree is  m = ';
  124. Mx3 = [Mx2,num2str(m)];
  125. Mx4 = 'The coefficients a(n) for cos(nx) are:';
  126. Mx5 = 'The coefficients b(n) for sin(nx) are:';
  127. clc,echo off,diary output,...
  128. disp(Mx1),disp(''),disp(Mx3),disp(''),...
  129. disp(Mx4),disp(''),for i=1:5:m+1,disp(A([i:min(i+4,m+1)]));end,...
  130. disp(Mx5),disp(''),for i=1:5:m+1,disp(B([i:min(i+4,m+1)]));end,...
  131. diary off,echo on
  132.  
  133. pause % Press any key to continue.
  134.  
  135. format short;
  136. points = [X;X;tp(A,B,X,m);Y-tp(A,B,X,m)]';
  137. Mx6 = '    x(k)      y(k)      P(x(k))   error';
  138. clc,echo off,diary output,...
  139. disp(''),disp(Mx6),disp(points),diary off,echo on
  140.  
  141.